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Pi  NONCONFORMING  FINITE  ELEMENT  METHOD  FOR  THE  SOLUTION  OF 
RADIATION  TRANSPORT  PROBLEMS 

KAB  SEOK  KANG* 

Abstract.  The  simulation  of  radiation  transport  in  the  optically  thick  flux-limited  diffusion  regime  has 
been  identified  as  one  of  the  most  time-consuming  tasks  within  large  simulation  codes.  Due  to  multimaterial 
complex  geometry,  the  radiation  transport  system  must  often  be  solved  on  unstructured  grids.  In  this  paper, 
we  investigate  the  behavior  and  the  benefits  of  the  unstructured  Pi  nonconforming  finite  element  method, 
which  has  proven  to  be  flexible  and  effective  on  related  transport  problems,  in  solving  unsteady  implicit 
nonlinear  radiation  diffusion  problems  using  Newton  and  Picard  linearization  methods. 

Key  words,  nonconforming  finite  elements,  radiation  transport,  inexact  Newton  linearization,  multigrid 
preconditioning 

Subject  classification.  Applied  and  Numerical  Mathematics 

1.  Introduction.  Radiation  transport  in  astrophysical  phenomena  and  inertial  confinement  fusion  is 
often  modeled  using  a  diffusion  approximation  [12,  17,  18,  20,  21,  22,  24].  When  the  radiation  field  is  not  in 
thermodynamic  equilibrium  with  the  material  a  coupled  set  of  time  dependent  diffusion  equations  is  used  to 
describe  energy  transport.  These  equations  are  highly  nonlinear  and  exhibit  multiple  time  and  space  scales. 
Implicit  integration  methods  are  desired  to  overcome  time  step  restrictions. 

Nonconforming  finite-element  methods  have  proven  flexible  and  effective  on  incompressible  fluid  flow 
problems  such  as  incompressible  Stokes  and  Navier-Stokes  equations  [10,  11].  In  the  Pi  nonconforming 
method,  the  degrees  of  freedom  lie  on  midpoints  of  edges.  Therefore,  the  number  of  connections  of  degrees 
of  freedom  with  each  others  at  most  four  (four  at  interior  edges  and  two  at  boundary  edges)  which  is  the 
same  number  of  connections  of  degrees  of  freedom  in  structured  finite  difference  methods.  In  contrast, 
in  the  Pi  conforming  method,  the  number  of  connections  of  degrees  of  freedom  is  at  least  four  except  at 
boundary  points,  and  depends  the  triangulation  and  position  of  points.  The  number  of  connections  of 
degrees  of  freedom  determines  the  number  of  nonzero  entries  of  generated  matrices  and  plays  an  essential 
role  in  performance  of  parallel  implementations  because  of  the  communication  required  in  kernel  operations 
like  matric- vector  multiplication.  Pi  nonconforming  methods  generate  matrices  that  have  a  constant  small 
number  of  nonzero  entries  for  each  row,  and  therefore  have  some  advantages  in  parallel  implementation  and 
performance. 

Because  many  nonlinear  elliptic  problems  are  well  solved  by  conforming  finite  element  methods,  non- 
conforming  methods  are  still  rare  for  such  problems.  However  nonconforming  methods  may  resolve  features 
of  solutions  of  nonlinear  problems  not  well  represented  by  conforming  methods.  In  this  research,  a  noncon¬ 
forming  methodis  shown  to  resolve  very  sharp  changes  of  energies  on  heterogeneous  domains.  The  results 
axe  very  similar  to  the  solutions  of  the  finite  volume  method  with  an  edge-based  flux  limiter  [19]. 

To  solve  nonlinear  problems,  one  usually  employs  linearization  techniques.  Many  modelers  use  Picard 
and  Newton  methods  to  linearize.  Picard’s  method  is  easy  to  understand  and  implement,  but  converges 
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also  partially  supported  by  postdoctoral  fellowships  program  from  Korea  Science  &  Engineering  Foundation  (KOSEF). 
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slowly.  Newton’s  method  has  a  second-order  convergence  rate  but  requires  the  Jacobian  of  the  original 
nonlinear  system.  In  many  nonlinear  problems,  an  inexact  Newton  method  works  well,  with  less  storage 
and  operation  count  expense  [8].  In  this  paper,  we  study  the  behavior  of  these  three  methods  on  a  model 
radiation  transport  problem. 

Because  the  system  generated  from  some  linearization  of  the  nonlinear  problem  is  usually  nonsymmetric, 
we  use  preconditioned  GMRES  [23].  As  a  preconditioner,  we  consider  multigrid.  Multigrid  represents  an 
important  advance  in  algorithmic  efficiency  for  the  solution  of  large  problems  [2,  3,  4,  14,  19,  25]. 

To  use  multigrid,  we  need  to  define  intergrid  transfer  operators  between  nonconforming  finite-element 
spaces.  Due  to  the  non-nestedness  of  nonconforming  spaces,  there  is  no  natural  intergrid  transfer  operator.  In 
previous  studies  of  the  nonconforming  multigrid  method,  the  average  value  of  two  adjacent  elements  is  used  to 
get  the  interpolated  value  at  a  node.  Nonconforming  multigrid  with  this  intergrid  transfer  operator  is  a  good 
solver  for  linear  systems  and  some  nonlinear  systems  with  smooth  nonlinear  coefficients  [1,  5,  6,  9,  15,  16]. 
However  this  intergrid  transfer  operator  does  not  preserve  positivity  of  functions,  which  is  an  essential  part  of 
radiation  transport  problems  because  energy  and  temperature  are  always  positive.  Therefore  some  nonlinear 
problems  with  discontinuous  coefficients,  bound  constraints  on  solutions,  and  rapidly  changing  solutions, 
like  the  radiation  transport  problem,  cannot  use  this  intergrid  transfer  operator  because  the  coarse  level 
approximation  obtained  from  the  fine  level  approximation  does  not  satisfy  solution  bounds,  and  one  cannot 
generate  the  coarse  level  systems  or  solve  the  coarse  level  problems  [15].  To  overcome  these  difficulties,  we 
use  a  new  and  simple  intergrid  transfer  operator  that  preserves  positivity  and  solves  the  above  mentioned 
problem.  However  multigrid  with  this  intergrid  transfer  operator  is  slower  than  with  the  previous  operator. 
Therefore,  we  use  the  simple  intergrid  transfer  operator  to  derive  coarse  level  systems  and  the  average  value 
intergrid  operator  to  solve  the  linear  systems. 

The  rest  of  the  paper  is  organized  as  follows.  In  section  2,  we  describe  a  model  radiation  transport 
and  its  Pi  nonconforming  discretization.  In  section  3,  we  consider  a  discretization  in  time,  derive  the 
linearizations  by  Picard  and  Newton  method,  and  describe  the  Inexact  Newton  method.  In  section  4,  we 
describe  preconditioned  GMRES  and  the  nonconforming  multigrid  preconditioner.  Numerical  experiments 
are  given  in  section  5. 

2.  Radiation  transport  model  and  Pi  nonconforming  discretization.  Under  the  assumption  of 
an  optically  thick  medium  (short  mean  free  path  of  photons)  a  first-principles  statement  of  radiation  transport 
reduces  to  the  radiation  diffusion  limit.  A  particular  idealized  dimensionless  form  of  such  a  system,  known 
as  the  “2T”  model,  can  be  written  as: 


with 


FiP 

—  -\/-(DrS/E)=aa{T^-E), 

(2.2.1) 

^-V-(AVT)  =  -(7a(T''-£), 

(2.2.2) 

(2.2.3) 

Here,  E{x^t)  represents  the  photon  energy,  T{x,t)  is  the  material  temperature,  da  is  the  opacity,  and 
K  is  the  material  conductivity.  In  the  non-equilibrium  case,  the  nonlinear  source  terms  on  the  right-hand 
side  are  nonzero  and  govern  the  transfer  of  energy  between  the  radiation  field  and  material  temperature. 
Additional  nonlinearities  are  generated  by  the  particular  form  of  the  diffusion  coefficients,  which  are  functions 
of  the  E  and  T  fields.  In  particular,  the  energy  diffusion  coefficient,  Dr(T,  E)  contains  the  term  |  VP|  which 
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Fig.  1.  Domain  of  model  problem 


refers  to  the  gradient  of  E.  This  limiter  term  is  an  artificial  means  of  ensuring  physically  meaningful  energy 
propagation  speeds  (i.e.  no  faster  than  the  speed  of  light).  The  atomic  number  z  is  a  material  coefficient, 
and  while  it  may  be  highly  variable,  it  is  only  a  function  of  position  (i.e.  z  =  f{x,y)  in  two  dimensions). 

The  two  model  problems  considered  in  this  study  are  taken  from  [19]  and  depicted  in  Figure  1.  We 
consider  a  unit  square  domain  of  similar  material  with  atomic  number  z  =  1  and  a  unit  square  domain  of 
two  dissimilar  materials,  where  the  outer  region  contains  material  with  an  atomic  number  of  ^  =  1  and  the 
inner  region  (1/3  <  a:  <  2/3,  1/3  <  y  <  2/3)  contains  material  with  an  atomic  number  oi  z  =  10.  The  top 
and  bottom  walls  are  insulated,  and  inlet  and  outlet  boundaries  are  specified  using  mixed  (Robin)  boundary 
conditions,  as  shown  in  the  figure.  For  convenience,  we  represent  the  boundary  a;  =  0,  0<2/<l  and  x  —  1, 
0  <  J/  <  1  by  To,  and  otherwise  by  Fi.  Then  the  boundary  condition  of  the  problem  is 


on  Fo, 
on  —  Fo, 
on  50, 

where  n  is  the  local  outward  normal  vector  of  the  boundary. 

Equations  (2.2.1)  and  (2.2.2)  form  a  system  of  coupled  nonlinear  partial  differential  equations  which 
must  be  discretized  in  space  and  time.  In  this  section,  we  consider  a  discretization  in  space  and  will  consider 
a  discretization  in  time  in  the  next  section. 

The  variational  form  of  (2.2.1)  and  (2.2.2)  can  be  written  as  follows:  Find  (jE?,T)  €  (iF^(O)nL^([0,T]))^ 
such  that 


E  1  dE  _ 

4  6(Ta  dn  ~ 

^  =  0 
dn  ’ 

^=0 
dn  ’ 


f  ^udx  +  f  DrVE  ■  Vudx  +  [  ^Euda 

Jq  Jq  JrQ  2 

—  I  aa{{T)^  —  E)udx  “  f  2gu(kr  =  0, 
Jq  Jtq 

f  ^vdx  +  f  DtVT  ■  Vvdx  +  f  aa{{Ty  -  E)vdx  =  0, 
Jq  Jq  Jq 


(2.2.4) 

(2.2.5) 


for  all  {u,v)  e  (fr^(O))^  and  for  all  t  e  [0,tmaxj* 

We  discretize  0  by  using  a  triangular  grid  containing  edges,  shown  in  Figure  2.  The  grid  is  generated  by 
connecting  of  the  midpoints  of  the  edges  of  the  triangles  from  the  coarsest  discretization  71,  which  contains 


Fig.  2.  Discretization  of  Domain 

edges  and  conforms  to  the  material  interface  boundaries  in  such  a  way  that  no  triangle  edges  cross  this 
boundary.  Let  hj  and  =  7j,  for  j  =  1, . . .  ,  J,  be  given,  where  7j  is  a  partition  of  Q  into  triangles  and  hj 
is  the  maximum  diameter  of  the  elements  of  Tj. 

Define  the  Pi -nonconforming  finite  element  spaces 

Vj  =  {v  €  :  v\k  is  linear  for  all  K  e  Tj, 

V  is  continuous  at  the  midpoints  of  interior  edges}. 

Then  the  nonconforming  finite  element  discretization  of  (2.2.4)  and  (2.2.5)  can  be  written  as  :  Find 
(Eh.Th)  G  (V>  X  [0,tmax])^  such  that 


^udx  +  ^  Dr{Ek,THWEh  ■  Vudx  + 

r 

[  ^EhUda 

JTo  ^ 
r 

(2.2.6) 

-  /  -  Ek)udx  - 

JQ 

/  2guda  =  0, 

Jto 

f  ^vdx+  [  Dt{En,Tk)vn  vvdx+  f  <Ta{n)iin)* 
Jq  ^  JQ  Jn 

“  Eh)vdx  =  0, 

(2.2.7) 

for  all  {u,v)  G  Vj  and  for  alH  €  [0,tmax]- 

In  above  equations,  to  perform  the  integration  in  space,  we  use  a  three-point  quadrature  rule  on  each 
triangle  in  7j.  Because  the  points  where  the  degrees  of  freedom  are  defined  and  the  quadrature  points  of 
triangle  are  the  same,  we  can  easily  compute  the  integration  on  each  triangle  and 

f  D{x)(f>i^kdx  =  (2.2.8) 

Jk  ^ 

for  all  basis  functions  (j>i  of  V  and  K  £Tj.  Also,  because  Vu  is  a  piecewise  constant  on  each  triangle  in  Tj 
for  all  6  G  V",  we  compute  |Vn|  needed  in  Dr  exactly. 

3.  Time  integration  and  nonlinear  iteration.  In  this  section,  we  consider  a  discretization  in  time 
and  three  nonlinear  iterations,  i.e.,  Newton,  Picard  and  inexact  Newton  iteration. 
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The  time  derivatives  are  discretized  as  first-order  backward  differences,  with  lumping  of  the  mass  matrix, 
leading  to  an  implicit  scheme  which  requires  the  solution  of  a  nonlinear  problem  at  each  time  step.  This 
approach  is  first-order  accurate  in  time,  and  is  chosen  merely  for  convenience,  since  the  principal  objective 
is  the  study  of  the  solution  of  the  nonlinear  system.  Higher  order  temporal  discretizations  are  demonstrated 
to  be  worth  while  in  [18]. 

To  solve  the  nonlinear  problem  (2.2.6)  and  (2.2.7)  we  consider  the  Picard  linearization  method  and  the 
Newton  linearization  method.  In  both  methods,  we  need  to  solve  linear  systems  to  get  corrections  at  eaeh 
nonlinear  iteration  step. 

The  fully  implicit  Picard  linearization  method  separates  the  operators  into  linear  parts  and  nonlinear 
parts  and  all  nonlinear  parts  are  evaluated  at  the  previous  nonlinear  iteration  level,  k  —  1.  This  results  in 
the  following  system  of  equations: 


i 

L 


En,k 

At 

udx  + 

f  D'!;:*-'^VEl*  -Vudx 
Jq 

r 

'^uda 

(3.3.1) 

/ 

Jq 

fc-1^3yn,fc  _  El'‘‘)udx  - 

I  2gud(T 
Jto 

=  0, 

rpTi^k 

rpn—l 

■vdx 

At 

p 

Jq 

(3.3.2) 

^n,*-l((yn,A=-l)3rn,.  _ 

E^''°)vdx  -- 

=  0, 

for  all  (u,v)  e  Vj.  Because  (3.3.1)  and  (3.3.2)  are  linear  systems  in  (J5^’*,T^’*'),  we  can  easily  calculate 
their  Jacobian. 

To  get  the  corrections  {6E,  ST)  in  the  Picard  Method  at  level  k,  we  solve  the  following  linear  systems. 


f  ^udx+  f  D”’^-'^V5E-Vudx+  [  IdEuda 

Jn^t  Jo  Jto^ 

-  f  -  SE)udx  = 

Jq 

^-^vdx+  f  D7'’°~'^V8T  -Vvdx 

A*  Jo 

+  f  -  SE)vdx  = 

Jq 


L 


for  ail  {u,  v)  e  Vj  where 


(3.3.3) 


(3.3.4) 


(3.3.5) 


(3.3.6) 


For  the  fully  implicit  Newton  linearization  method  it  is  somewhat  more  complicated  to  compute  the 
Jacobian  at  approximate  solution  points.  To  get  the  Jacobian,  we  have  to  calculate  the  derivatives  of  the 
system  with  respect  to  (0i,^t)  for  all  basis  functions  in  Vj  xVj. 
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As  the  result  of  differentiation  with  respect  to  to  get  the  corrections  {5E,ST)  in  Newton’s 

Method  at  level  A;,  we  solve  the  following  linear  systems. 

SE 


.  S/udx 


f  ^udx+  I  D”''‘~^VSE  Vudx+  [ 

Jii  Ju  Jq 

+  /  D^’^-^STVE]l;’'‘-^  ■  Vudx  +  /  \5Euda 
Jq  ’  Jvo  2 

-  ^  ^1  +  3^^^  6Tudx  +  £  a2’'‘-^SEudx  = 


i 


5T  ^ 
—vdx  + 

Jn 


+ 


/  A"’*"' 

Jq 

/„<■“•  (' 


V^T  •  Vudar  + 


\n,k—l 


nTlyk  —  l 


STvdx  ■ 


Jq 

f 

Jq 


(SEwdx  =  F^''‘-\v), 


for  all  (u,  r)  €  where 


r.T  ’ 

Of  ^  pin  ,fc  I 

where  *  can  be  easily  evaluated  on  each  triangle  in  7j. 

o{Eh'  )i 

After  linearization,  we  have  to  solve  the  linear  systems 


jk—l 


SE 

ST 


-  rE 


rtTlyk  —  i 


(3.3.7) 


(3.3.8) 


(3.3.9) 


for  each  step  where  J*  ^  is  a  Jacobian,  which  is  computed  by  Picard’s  method  or  Newton’s  method. 
In  either  method,  we  need  for  robustness  to  control  the  step  length  a  where 


In  this  study,  we  control  the  step  length  by  simply  halving  a  until  the  residual  of  the  updated  solution  is 
less  than  the  previous  residual.  In  this  control,  we  sometimes  fail  to  get  a  proper  step  length,  so  we  stop  at  a 
fixed  step  length  and  perform  the  next  nonlinear  iteration.  If  the  number  of  failures  exceeds  a  fixed  number, 
then  we  go  to  next  time  steps  by  using  the  best  approximation,  which  has  the  smallest  nonlinear  residual. 

Remark  3.1.  The  Newton  method  has,  asymptotically,  a  second  order  convergence  for  nonlinear  prob¬ 
lems  and  the  Picard  method  has  only  a  first  order  convergence.  However  the  resulting  linear  problem  of  the 
Picard  method  is  more  easily  solved  than  that  of  the  Newton  method  because  the  Picard  method  lacks  the 
convection  term  as  described  in  ref  [7]. 

To  improve  the  efficiency  of  the  Newton  method,  we  can  use  an  inexact  Newton  method  [8].  When 
the  Newton  iteration  is  “far”  from  convergence  (i.e.,  the  residual  is  large)  there  is  no  reason  to  solve  the 
linear  system  accurately.  However,  when  the  Newton  iteration  is  “close”  (i.e.,  the  residual  is  “small”)  the 
convergence  rate  of  Newton’s  method  is  tightly  coupled  to  the  accuracy  of  the  linear  solution.  To  adjust  the 
amount  of  work  done  in  the  linear  solve  (via  a  convergence  tolerance)  we  employ  an  inexact  Newton  method. 


6 


In  the  inexact  Newton  approach,  the  convergence  criteria  for  the  linear  solver  is  proportional  to  the  residual 
in  the  nonlinear  iteration.  In  equation  form  this  is 


J 


fe-i 


(3.3.10) 


where  72  =  1-0  x  10”^  is  the  value  used  in  this  study  unless  otherwise  noted.  We  note  that  [13]  shows  how 
to  adaptively  select  72  to  recover  asymptotically  full  second  order  convergence. 


4.  PGMRES  and  multigrid  preconditioning.  In  this  section,  we  explain  PGMRES,  which  is  a 
combination  of  a  Krylov-based  linear  iterative  method,  and  multigrid,  which  is  well  known  as  a  successful 
preconditioner,  as  well  as  a  scalable  solver  even  in  unaccelerated  form,  for  many  problems. 

GMRES  [23]  is  a  well  known  solver  for  non-Hermitian  problems.  In  practice,  GMRES  can  be  restarted 
after  m  steps,  where  m  is  some  fixed  integer  parameter,  to  save  storage  by  accepting  a  generally  less  rapid 
convergence. 

We  describe  the  restarted  PGMRES  for  solving 


Ajx  —  h 


(4.4.1) 


with  preconditioning  matrix  Bj. 

PGMRES(m)  Algorithm  4.1. 

(1)  Start  :  Choose  xq  and  compute  ro  =  Bj{b  -  AjXq),  /?  =  Urolb  and  ui  =  rof/S. 

(2)  Iterate  :  For  j  =  1, . . .  ,  m  do: 

Compute  w  :=  BjAjVj 

For  i  =  1,...  , j,  do: 

hi,j  ■■=  (w,Vi) 

w:=w  -  jijVi 

Enddo 

Compute  hj^i^j  =  ||«;||2  and  Uj+i  =  w/hj^ij 
Enddo 

(3)  Form  the  approximation  solution: 

Define  Vm  :=  [vi,..-  ,Vm], 

and  set  Xm-Xo^  Vmym,  where  ym  minimizes  ||/?ei  -  Hmy\l  y  ^ 

(4)  Restart: 

Compute  Tm  =  Bj(6  -  AjXm);  if  satisfied  then  stop 

else  compute  xq  :=  Xm,  P  —  Ikmll  and  vi  =  TmlP  and  go  to  (2). 

Arnold!  iteration  constructs  an  orthogonal  basis  of  the  left  preconditioned  Krylov  subspace 

Span{ro,  BjAjv^^, . . .  ,  (BjAj)^"^ro}. 

It  uses  a  modified  Gram-Schmidt  process,  in  which  the  new  vector  to  be  orthogonalized  is  obtained  from 
the  previous  vector  in  the  process.  All  residual  vectors  and  their  norms  that  are  computed  by  the  algorithm 
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correspond  to  the  preconditioned  residuals,  namely,  Zm  =  Bj{h-  AjXm)^  instead  of  the  original  (unprecom 
ditioned)  residual  b  —  AjXm-  In  addition,  there  is  no  easy  access  to  these  unpreconditioned  residuals,  unless 
they  are  computed  explicitly.  So  we  monitor  these  preconditioned  residuals  to  stop  PGMRES  iteration  to 
solve  linear  problem. 

Next,  we  consider  Multigrid  Preconditioner  Bj, 

To  define  a  multigrid  method,  we  need  to  define  intergrid  transfer  operators  between  nonconforming 
finite  element  spaces.  Due  to  the  non-nestedness  of  nonconforming  spaces,  there  is  not  a  natural  intergrid 
transfer  operator.  In  previous  studies  of  nonconforming  multigrid  method[l,  5,  6,  9],  average  value  of  two 
adjacent  elements  are  used  to  set  the  value  of  a  node.  A  nonconforming  multigrid  method  with  this  intergrid 
transfer  operator  is  a  good  solver  for  linear  systems  and  some  nonlinear  systems  that  have  smooth  nonlinear 
coefficients. 

To  get  the  coarse  level  approximate  linear  system  for  (3.3.9),  we  need  coarse  level  approximations  of 

jf  approximate  solution  varies  rapidly  in  space, 

thensomecoarselevelapproximationsof(E^’^”^,T^’^“^)  may  have  negative  values.  However 
axe  required  to  be  positive  for  the  computation  of  Either  we  cannot  generate  the  coarse  level  systems 

or  they  may  become  nearly  singular,  making  it  hard  to  solve  the  coarse  level  problems. 

To  overcome  these  difficulties,  we  use  a  new  and  simple  intergrid  transfer  operator  called  the  covolume- 
based  intergrid  transfer  operator,  which  preserves  only  piecewise  constant  functions  [15].  It  is  well  known 
that,  to  get  a  good  convergence  factor  in  multigrid  algorithms,  intergrid  transfer  operators  should  preserve 
higher  order  functions  [19].  Therefore  the  multigrid  method  with  this  intergrid  transfer  operator  converges 
slowly  compared  to  average  value  intergrid  operator  to  solve  linear  systems.  However  preservation  of  pos¬ 
itivity  of  nodal  values  of  the  fields  is  critical.  So,  we  use  the  covolume-based  intergrid  transfer  operator 
to  obatain  the  coarse  level  systems  and  the  average  value  intergrid  operator  to  interpolate  the  solution  be¬ 
tween  levels  (coarse-to-fine  and  fine-to-coarse)  when  solving  the  linear  systems  in  Picard’method  or  Newton’s 
method. 

Let  Aj\  (Vj)^  j  =  1, . . .  ,  J  be  the  discretization  operator  on  level  j  and  Ij  : 

j  =  2, . . .  ,  J,  be  the  coarse-to-fine  intergrid  transfer  operator.  Also,  we  define  the  fine-to-coarse  intergrid 
transfer  operator  PjLi  : 

{IjV,w)  =  Vv  e  €  {Vjf. 

Finally,  let  Rj  :  (Vj)^  -*■  (Vj)^  for  j  =  1,...  ,  J  be  the  linear  smoothing  operators,  let  Rj  denote  the 
adjoint  of  Ej  with  respect  to  the  (•,  •)  inner  product,  and  define 

I  odd, 

^  ^  ®ven. 

Following  [2],  the  multigrid  operator  Bj  :  (Vj)^  (^j)^  is  defined  recursively  as  follows. 

Multigrid  Algorithm  4.2.  Let  1  <  j  <  J  and  p  be  a  positive  integer.  Set  Bi  =  A^^.  Assume  that 
Bj^i  has  been  defined  and  define  Bjg  for  p  €  (Vj)^  by 

(1)  Set  X®  =  0  and  —  0. 

(2)  Define  for  /  =  1, . . .  ,m(j)  by 


8 


(3)  Define  4-  where  for  i  =  1, . . .  ,p  is  defined  by 

q*  =  q'-^  +  B^_i[P9_i(5  -  -  Aj-iq^-^]. 

(4)  Define  for  I  =  m{j)  +  1, . . .  ,  2m(j)  by 

y^  =  y^~^  +  -  Ajy^~'^), 

(5)  Set  Bjg  =  y^^(^\ 

In  Multigrid  algorithm  4.2,  m{j)  gives  the  number  of  pre-  and  post-smoothing  iterations  and  can  vary 
as  a  function  of  j.  If  p  =  1,  we  have  a  F-cycle  multigrid  algorithm.  If  p  =  2,  we  have  a  VF-cycle  multigrid 
algorithm.  Other  versions  of  multigrid  algorithms  without  pre-  or  post-smoothing  iterative  can  be  analyzed 
similarly.  A  variable  F-cycle  multigrid  algorithm  is  that  for  which  the  number  of  smoothing  m(j)  increases 
exponentially  as  j  decreases  (i.e.,  p  =  1  and  7n(j)  —  2'^“^). 

Remark  4.1.  One  can  use  the  multigrid  algorithm  to  solve  the  systems  as  a  free-standing  iterative 
method.  Usually,  one  uses  V -cycle  and  W -cycle  multigrid  algorithms  to  this  end  and  uses  V -cycle  and 
variable  V -cycle  multigrid  method  as  preconditioners  of  Krylov-type  methods  such  as  PCG,  because,  when  Aj 
is  symmetric  positive  definite,  the  V -cycle  multigrid  operator  Bj  is  a  symmetric  positive  definite  operator  on 
{Vj)^,  but  the  W -cycle  multigrid  operator  is  not  in  generally  [3].  Many  researchers  show  that  convergence 
of  W -cycle  multigrid  for  the  nonconforming  and  conforming  cases  and  V -cycle  multigrid  for  the  conforming 
case  are  good  preconditioners  [1,  2,  5,  6,  9,  14,  16,  25].  In  this  problem,  we  use  V -cycle  multigrid  method 
as  a  preconditioner  of  GMRES. 

5.  Algorithm  performance  and  results.  In  this  section,  we  study  the  performance  of  the  Newton, 
Picard,  and  inexact  Newton  methods  on  Pi  nonconforming  finite  element  method  on  two  model  problems 
with  the  only  difference  between  the  problems  being  homogeneity.  In  the  two  examples,  we  use  the  same 
triangulations,  namely  12800  triangles,  19296  edges,  and  6497  vertices.  Because  nodes  are  on  midpoints  of 
edges  in  a  Pi  nonconforming  method,  the  number  of  degrees  of  freedom  of  this  problem  is  38592. 

For  problem  1,  we  consider  a  homogeneous  material  with  atomic  number  z  —  I  and  k  =  0.01  on  the 
whole  domain.  The  initial  conditions  are  P®  =  1.0  x  10”®  and  T®  =  The  problem  is  run  out  to 

time  t  =  3.0  and  nonlinear  convergence  tolerance  within  a  time  step  is  defined  as  ||F(u*^)|l2  <  1.0  x  10”®  for 
problem  1.  We  run  with  several  time  steps  of  0.001,  0.002,  0.005,  and  0.01. 

In  Figure  3,  we  plot  the  contours  of  temperature  T  at  t  —  1.0, 2.0, 3.0.  Table  1  compares  linear  solve 
requirements  and  nonlinear  iterations.  Figure  4  depicts  the  nonlinear  convergence  behavior  of  Newton 
method,  Picard  method,  and  Inexact  Newton  method  at  time  t  =  1.0. 

Figure  3  shows  that  contours  of  temperature  propagate  parallel  to  the  inlet  boundary  and  reproduce  on 
an  unstructured  grid  the  propagation  of  the  one-dimensional  case.  Table  1  and  Figure  4  show  that  Newton’s 
Method  is  very  efficient  compared  to  Picard’s  method,  and  slightly  more  efficient  compared  to  the  inexact 
Newton  method,  in  terms  of  nonlinear  iterations  per  time  step. 

Inexact  Newton  needs  more  nonlinear  iterations  in  comparison  to  Newton’s  method,  but  has  the  best 
performance  overall  because  this  method  needs  the  smallest  number  of  linear  iterations.  Also,  Table  1  shows 
that  the  number  of  linear  iterations  in  each  nonlinear  iteration  of  the  Picard  method  is  smaller  than  that  of 
the  Newton  Method.  This  means  that  the  linear  systems  from  Picard’s  method  are  more  easily  solved  than 
the  linear  systems  from  Newton  method. 

In  Table  2,  we  report  the  accuracy  as  a  function  of  time  steps  by  the  -error  of  the  solution  which 
is  defined  as  ||it  —  «base||2  where  Ubase  is  obtained  by  using  a  time  step  0.0001.  This  result  shows  that  the 
L^-error  in  time  is  first  order. 


9 


(a)  t  =  1.0 


(b)  t  =  2.0 


(c)  t  —  3.0 

Fig.  3.  Contour  of  Temperature  of  Problem  1 

For  problem  2,  we  consider  an  inhomogeneous  material  with  atomic  number  z  =  10  inside  the  box  and 
.2:  =  1.0  outside,  as  shown  in  Figure  5.  We  changed  the  nonlinear  convergence  tolerance  within  a  time  step 
to  be  ||F(u^)||2  <  1.0  x  10“^  to  reduce  the  simulation  times. 

In  Figure  6,  we  plot  the  contour  of  temperature  T  at  f  =  1.0, 2.0, 3.0, 4.0, 5.0.  Table  3  compares  linear 
solve  requirements,  nonlinear  iterations,  and  number  of  failures  to  meet  the  convergence  tolerance.  Figures  8, 
9,  10  demonstrate  the  nonlinear  convergence  behavior  of  the  Newton,  Picard,  and  inexact  Newton  methods 
at  times  t  =  1.0,  t  =  2.5,  and  t  =  4.0. 

As  energy  propagates,  temperatures  rapidly  change  near  the  front  and  near  the  layer  where  the  two 
different  materials  meet.  As  more  time  passes,  the  temperature  smoothly  propagates.  Figures  8,  9, 10  show 
that  there  are  many  step  length  controls  to  get  the  solution  of  the  nonlinear  problem  when  the  solution 
changes  rapidly  {t  =  1.0,  t  =  2.5)  but  there  is  no  need  for  step  length  control  when  the  solution  is  smooth 
(t  =  4.0)  in  any  of  the  three  methods. 
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Table  1.  Algorithm  performance  as  a  function  of  time  step  for  problem  1  and  time  period  of  3.0 


Method 

dt 

of  dt 

tot  # 

nonlin¬ 

ear 

ave  # 

nonlin 

iter  /dt 

tot  # 

linear 

ave  # 
lin  /dt 

ave  #  lin 
/nonlin 

0.001 

6226 

2.1 

49851 

16.6 

Newton 

0.002 

4116 

2.7 

38970 

26.0 

9.5 

method 

0.005 

600 

2120 

3.5 

28197 

46.8 

13.3 

0.01 

300 

1334 

4.4 

21986 

73.7 

16.5 

0.001 

3000 

24935 

8.3 

181227 

60.4 

7.3 

Picard 

0.002 

1500 

15389 

10.3 

126152 

84.1 

8.2 

Method 

0.005 

600 

7784 

13.0 

70165 

116.9 

9.0 

0.01 

300 

5320 

17.7 

46194 

154.0 

8.7 

Inexact 

0.001 

3000 

8648 

2.9 

28468 

9.5 

3.3 

Newton 

0.002 

1500 

4534 

3.0 

15928 

10.6 

3.5 

Method 

0.005 

600 

2254 

3.8 

9878 

16.6 

4.4 

0.01 

300 

1450 

4.8 

7761 

25.9 

5.4 

Table  2.  -error  att  —  3.0 


time  steps 

(error) 

0.001 

0.00884 

0.002 

0.01796 

0.005 

0.04060 

0.01 

0.06676 

Figure  6  shows  that  the  solution  of  the  nonconforming  finite  element  method  is  very  similar  to  the 
solution  of  finite  volume  method  with  edge-based  flux  limiter  [20]. 

In  the  aspect  of  performance,  the  behavior  of  problem  2  is  similar  to  problem  1  with  the  exception  that 
problem  1  does  not  require  step  length  control. 

To  estimate  the  accuracy  as  a  function  of  time  step  size,  we  report  the  L^-error  of  the  solution  in  Table 
4  (based  on  an  accurate  solution  with  dt  —  0.0001).  The  relative  L^-error  of  simulations  with  dt  =  0.002, 
0.005,  0.01  compared  to  the  L^-error  of  dt  =  0.001  is  ploted  as  a  function  of  time  in  Figure  7.  These  results 
show  that  the  L^-error  in  time  is  first  order  at  the  beginning  of  simulation  until  t  =  3.0  but  gradually 
deteriorates.  This  deterioration  may  be  introduced  by  the  nonlinear  convergence  error  within  a  time  step 
because  the  accumulation  of  the  nonlinear  convergence  error  will  dominate  other  errors  (space  and  time 
discretization  error)  as  time  steps  grows.  If  we  use  a  finer  nonlinear  convergence  tolerance,  then  we  can 
delay  this  deterioration  to  longer  time. 

6.  Conclusions.  We  solved  unsteady  implicit  nonlinear  radiation  diffusion  problems  by  an  unstruc¬ 
tured  Pi  nonconforming  finite  element  method.  Pi  nonconforming  finite  element  methods  resolve  very  sharp 
changes  of  energies  on  the  heterogeneous  domains,  similarly  to  results  of  the  finite  volume  method  with  an 
edge-based  flux  limiter.  The  inexact  Newton  method  has  the  best  performance  overall  and  Preconditioned 
GMRES  with  nonconforming  multigrid  preconditioner  to  solve  linear  problems  works  well.  In  Pi  noncon¬ 
forming  multigrid,  the  covolume-based  intergrid  transfer  operators  are  useful  to  solve  radiation  transport 
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(a) 


Fig.  4.  Convergence  plot  on  problem  1  at  time  t  =  1.0,  (a)  Newton  Method,  (h)  Picard  Method,  (c)  Inexact  Newton  Method 


problems  because  the  positivity  preserving  property  is  needed. 
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Table  3.  Algorithm  performance  as  a  function  of  time  step  for  problem  2  and  time  period  of  5.0 


2 


i ^ _ , _ j 

4  6  8  10  12  iterations 

(a)  dt  =  0.001 


(c)  dt  =  0.005 


—  Newton  Method 
Picard  Method 

—  Inexact  Newton  Method 


Fig.  8.  Convergence  plot  at  time  t  =  1.0 
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—  Newton  Method 
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Fig.  10.  Convergence  plot  at  time  t  =  4.0 
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